###Import emissions parameters
 
import XLSX
df = XLSX.readxlsx("aux_params/parameters.xlsx")["emissions_total"]
E = df[:][2:end,2:3]

###Set input and output folder

if conduct=="pc"
    inputs_folder=string("results_pc")
    results_folder=string("results_pc/mps")
else
    inputs_folder=string("results_ic")
    results_folder=string("results_ic/mps")
end

###Compute MPS

QU=readdlm(string(inputs_folder,"/ur/Q.csv"),',',Float64)
QT=readdlm(string(inputs_folder,"/dt/Q.csv"),',',Float64)

mu_E = mean(E, dims=1)
sigma_E = std(E, dims=1)

b_grid = round.(collect(0:0.2:3), digits=1)
n_b = size(b_grid,1)
results = zeros(n_b,2)
counter = collect(1:n_b)

for n in counter

    b = b_grid[n]
    a_shift = (1-b).*mu_E
    Em = a_shift.+ b.*E

    mu_Em = mean(Em, dims=1)
    sigma_Em = std(Em, dims=1)

    dsigma = b^2
    dE_diff = (sum((QT-QU).*Em))
    
    results[n,1]= dsigma
    results[n,2]= dE_diff

end

writedlm(string(results_folder,"/dE.csv"), results,',')